Hyperdynamics for entropic systems: time-space compression and pair correlation 

function approximation 
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We develop a generalized hyperdynamics method, which is able to simulate slow dynamics in 
atomistic general (both energy and entropy-dominated) systems. We show that a few functionals of 
the pair correlation function, involving two-body entropy, form a low-dimensional collective space, 
which is a good approximation that is able to distinguish stable and transitional conformations. A 
bias potential, which raises the energy in stable regions, is constructed on the fly. We examine the 
slowly nucleation processes of a Lennard- Jones gas and show that our new method can generate 
correct long time dynamics without a prior knowledge. 
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Molecular dynamics (MD) simulations are typically 
limited to a time scale of less than a microsecond, so 
many interesting slow processes in chemistry, physics, bi- 
ology and materials science cannot be simulated directly. 
Recently, new methods, including kinetic Monte Carlo, 
transition path ensemble methods, minimal action/time 
methods and the constrained MD simulation 0, 0, 0, , 
have been developed to study the slow processes (for a 
review see [fj ) . They all require a prior knowledge of the 
system, which is often hard to obtain, and they can only 
deal with a few special processes inside a small part of 
the configurational space of the system. 

For many systems, the interesting dynamics are gov- 
erned by the infrequent, fast transitions between meta- 
stable regions; yet the systems spend most of their time 
in the stable regions, whose dynamics can be well de- 
scribed by some time-averaged properties. Hence we 
would coarse grain the stable configurations while keep- 
ing the needed details in the unstable regions that de- 
fine the transitions. Hyperdynamics, as developed by 
Voter 6] is an example of such a coarse-graining method. 
The hyperdynamics method treats the potential wells as 
the stable conformational regions that are separated by 
the saddle regions. A bias potential is designed to lift 
the energy of the system in these wells, while keeping 
the saddle regions intact. Dynamics on the biased po- 
tential leads to accelerated evolution from one stable re- 
gion to another. Based on transition state theory, the 
realistic escape time t rea i from the wells can be repro- 
duced, t rea i = At J2 i exp[/3AV(r(i;))], where At is the 
time step of MD, AV(r) is the applied bias potential 
along the simulated trajectory r(t). This method has 
been applied successfully to systems in which the rel- 
evant states correspond to deep wells in the potential 
energy, with dividing surfaces at the energy ridge tops 
separating these states ■ It has however not been clear 
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how to apply hyperdynamics in cases where the transi- 
tions are dominated by entropic considerations. In such 
cases, the potential energy alone is not enough to dis- 
tinguish stable and transition regions since some confor- 
mations with similar energy might belong to stable and 
transition regions, respectively. A complication that oc- 
curs even when trying to apply hyperdynamics to solids 
with fairly clearly defined stable regions is that after ap- 
plying the bias potential, the energy landscape becomes 
much flatter and the system can start to have entropic- 
like characteristics. These effects limit the improvement 
in the simulation rate that can be achieved by the hyper- 
dynamics method over the direct MD approach. Thus, 
although there are some attempts in the literature [H, 
to apply hyperdynamics to enhance conformational sam- 
pling in bio-systems, generally, accurate slow dynamics 
or kinetic rates can only be expected for relatively simple 
solids or low dimensional systems. 

In this Letter, we derive a more general hyperdynamics 
method that can be used to access longer time scales in 
fluids. We present explicit conditions for applying this 
method. We use the pair correlation function as a reliable 
means of identifying the important conformations and 
then construct the appropriate bias potentials in a lower- 
dimensional space while the simulation runs. We then 
examine the performance of this method by looking at 
the slow gas-liquid transition in a system of N identical 
Lennard-Jones particles. 

We begin this process by introducing a time- 
compressing transformation, dr — a(r)dt, where dr is 
a pseudo-time step, dt is the real time step, and the local 
dimensionless compression factor is given by a confor- 
mational function a(r), which is < 1. Thus the trajec- 
tory, r(t), can be rewritten as r(r) = r(r(i)) in a shorter 
pseudo-time interval, r = t J dr D(r;r(t);t) a(r), where 
D(r\ r(t); t) is the distribution of r(t) in the interval [0, t]. 
The compressed trajectory r(r) satisfies a new equation 
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of motion, 

— P =-— + V^P •— InoM (1) 

where P^ = MdRi/dr, M is the transformed mass, equal 
to ma 2 (r), and the m is the real mass of the particles. 
V(r) is the potential energy, and r is the simple dona- 
tion of the position vectors R; (i = 1, • • -N) of all the N 
particles. At first glance it would not appear to be ad- 
vantageous to directly generate r(r) from Eq.Q as very 
short time steps are necessary due to the small values of 
M . However, if we only focus on the long-time dynam- 
ics, we can use a smoother pseudo-trajectory 1Z(t) to 
replace r(r), provided that one requires the reproduced 
time from 1Z(t) to be the same as that from r(r). Thus, 
a sufficient condition to replace r(r) with 1Z(t) is that 
their distributions are the same. 

In general, the distribution along a finite-length trajec- 
tory is not easily known. However, in the time-consuming 
regions (stable regions), similar conformations would be 
visited many times even during a finite simulation time. 
Thus we can assume the distribution can be approxi- 
mated by D(r;r(t);t) oc exp(— /3V(r)). Many methods 
might be used to generate 1Z(t) with the required distri- 
bution. One simple method is to use a realistic trajectory 
corresponding to the local equilibrium of a new potential 
U(r) = V(r) — ksT In a(r). Actually, based on the same 
local equilibrium supposition, if one replaces the kinetic 
energy term of Eq.Q by its ensemble averaged value, 
< PjPj/M >= ksTSij, Eq.QJ is indeed the equation 
of motion of the particles with smaller mass M under 
the new potential U(r). Outside of the time-consuming 
regions, we choose a(r) as unity so that the MD time 
is realistic. If we select small a(r) only in the poten- 
tial well regions, we effectively have the hyperdynamics 
presented by Voter [(| . The condition for the use of com- 
pressed time is that the simulated trajectory under U (r) 
has the equilibrium distribution in all the biased (time- 
compressed) regions. Simply, we can suppose that the 
distribution is locally given by the Boltzmann one in the 
conformational regions where the value of the distribu- 
tion is larger than a critical value. In comparison with 
the transition state theory in original hyperdynamics im- 
plementation, the new approach makes it is easier to both 
determine the proper regions to bias and to design the 
appropriate bias potential. 

In solids, as the number of conformations inside sta- 
ble regions (potential energy wells) is often small and 
the distribution of long-time trajectories can reach lo- 
cal equilibrium, the hyperdynamics method works very 
well. However, in entropy-dominated systems (e.g. gas 
and liquids), the number of conformations (entropy) in- 
side time-consuming regions may be huge, and thus the 
Boltzmann distribution may not be reached in finite MD 
time. In this case, we average the neighboring distri- 
bution of the trajectory, D(r; a) — J D(r)dr, where 
a is the size of selected neighbors. For the smoothed 
distribution, in time-consuming regions, we can expect, 



D(r;a) oc f exp(— (3V(r))dr. Thus, by selecting smooth 
a(r) functions, the realistic time propagator can still be 
reproduced from the physical trajectory of U(r), pro- 
vided the averaged distribution of the U(r) trajectory 
satisfies the equation. 

As an illustration, we compress the 3A-dimensional 
flat conformational space r to a curved 3A-dimension 
q space, dq = A(r)dr, where A(r) = ^ is the Jaco- 
bian of the transformation. In q space, the trajectory 
q(t) = q(r(t)) satisfies a new equation of motion with 
a positive symmetric mass matrix M = B T (r) m B(r), 
where B T denotes the transpose, m is the real mass di- 
agonal matrix, and B = A -1 . Not losing any generality, 
we choose a diagonal M by rotating the dq space, and 
the new equation of motion is, 

d dV v _g WnMM m 

dt dq Y 2M J d i 

where P and Pj are the generalized momentum and its 
j component, respectively, and Mj is the j diagonal ele- 
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ment of M. If one replaces jj- with /c^T, Eq.@ is the 
equation of motion of heavier particles (mass Mj) un- 
der a new potential, W(q) — V(q) + fc^Tlog J(q), where 
J(q) = det(A) is the determinant of the Jacobian A of 
the transformation. Here, W(q) is indeed the free energy 
of system in q space. Thus, by inhomogeneously com- 
pressing the conformational space, we transform the orig- 
inal entropy-dominated r space to a energy-dominated q 
space with effective potential W(q). If the trajectory 
q{t) of W(q) corresponds to the Boltzmann distribution 
in some time-consuming regions, the corresponding tra- 
jectory r(q(t)) represents a local equilibrium in r space. 
Thus we can bias the effective potential W(q) in q space 
to extend the MD time scale. Actually, this can be done 
directly in r space without using the explicit transforma- 
tion A(r). 

The keys to successfully apply hyperdynamics are dis- 
tinguishing the conformations and designing suitable bias 
potentials AV(r) for the entire conformation space r. 
Obviously, AV(r) should have the same symmetry as 
V(r). Considering a simple case of N identical particles, 
we rewrite the conformational vector {Rj} {j — 1, N) 
as a density field, p(x) = J2j$( x ~ ^-i)- Here both x 
and R are the normal 3-dimcnsion spatial vectors. Since 
the neighboring conformations are identical in the view- 
point of slow dynamics, p(x) can be averaged to get a 
smooth function, p(x), by for example using a Gaussian 
function to replace the Dirac-<5 function. If the width 
of the Gaussian function is small, p(x) can be used to 
identify different conformations. Here we used a func- 
tional space to replace the SA^-dimension conformation 
space, but actually, the physically allowed p(x) only oc- 
cupies a very small part of the functional space. By 
neglecting multi-body correlations and directional cor- 
relations, we can approximate the density field p{x) (or 
conformations) by using some bin-averaged values of the 
radial pair correlation function g{x) of the conforma- 
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tions (here, x donates the length of the spatial vector x), 
gi = x /a g( x + x i)dx. Here gi is an average value of g(x) 
in a small bin (x,— A/2, x,*+A/2). Thus, each conforma- 
tion corresponds to the group of gi that defines a point in 
5 space. The spatial neighbors of the conformation and 
their symmetric companions will also correspond to the 
same g point. If the bin size is very small, all conforma- 
tions with the same {gi} are identical in the slow dynam- 
ics viewpoint and thus the bias potential of hyperdynam- 
ics can be written as a function in the low-dimension g 
space, AF(R A ') = f({g x (R N )}). To better identify con- 
formations even when larger bin sizes are used, some im- 
portant dynamics-related physical variables, such as the 
potential energy ^(R^), can be added into the {gi} vari- 
able group. Another important variable is the two-body 
entropy, S2 — —2irp J[g(x) In g(x) — g(x) + l]x 2 dx. which 
forms the main part (about 90 percent) of the macro- 
scopic excess entropy ^lj- Similarly, it is also possible 
to use some other functional of g(x) to replace some gi. 
In special systems, it may be useful to add some spe- 
cial order parameters Oj to take into account possible 
multi-body correlations and thereby decrease the needed 
number of gi. Finally, we have a group (of order 10) of 
general collective variables denoted as 5 — {5 J }, which 
might involve V, S2, some {gi} and some possible {Oj}, 
to identify conformations and form an appropriate bias 
potential. In general, we construct the bias potential as 
AV(S(R N )) = k B Tf+(\nD(S)/D c ), where, D(S) is the 
distribution of the simulated trajectory, D c is the selected 
critical value. The function f+(z) = z for larger positive 
z, and smoothly approaches as z decreases to 0. The 
designed bias potential will generate a flatter distribu- 
tion in 5 space while the biased regions are still visited 
often enough for the system to reach local equilibrium. 
The bias potential can be formed gradually. First, we 
generate a long non-biased MD trajectory to calculate 
D(S) and form a (small) bias in some 5 regions. Next a 
long trajectory is simulated in the biased system, which 
in turn generates the basis for the next bias. This pro- 
cess is repeated until the desired long realistic time is 
reached. Thus, we can gradually study dynamics at ever 
longer time scales, but at the expense of the details seen 
in fast dynamics. With the biased potential, the bias 
force on particles can be calculated by the chain rule of 
differentiation, Af, = -£\ ^srjg. 

We have examined the general hyperdynamics method 
in a simple system of N identical Lennard- Jones (LJ) 
particles and studied the slow gas/liquid transition in the 
NVT ensemble. We used the truncated and shifted LJ 
potential with r c = 2.5, and the reduced units: ey = 1, 
cry = 1 and the mass of particles m — 1. The tempera- 
ture is fixed as T — 0.613. The velocity Verlet algorithm 
was used to integrate the Langevin equation of motion. 
Obviously, in such a system, the potential energy of tran- 
sitional conformations ( liquid drops with critical size) is 
lower than that of the stable gas phase. Thus, some sim- 
ple bias methods, for example, just lifting the energy of 
all lower-energy conformations [lOj , cannot work at all in 



such a system. Actually, for general entropy-important 
systems, using the potential alone is not enough to iden- 
tify the transitional conformations and then form the 
bias potential. Only in some special low-dimension sys- 
tems |9J where the entropy is occasionally not important, 
exceptions may be found. 

We initially studied a gas phase with a relatively large 
saturation ( a number density n — 0.02 for a system of 
N = 1000 particles) for which a liquid drop forms in nor- 
mal MD without any bias. The phase transition happens 
in a narrow time window where the potential V, the two- 
body entropy 1S2 and the pair correlation function g{x) 
are found to change drastically in the gas/liquid phase 
transition process. However, inside each small S% range, 
we found g(x) only fluctuates slightly around its average 
value. This shows that the two-body entropy £2 inte- 
grates the information of the g(x). Thus, in this simple 
system, we use only two functionals of g(x), namely V 
and 52, to form the collective space while designing the 
bias potential. By examining the <7,-, which are available 
from our simulations, we were able to show that V and 
52 can sufficiently identify transition states and stable 
states, at least in this simple system. Since 52 gives the 
main part of the entropy, we expect that 52 and V are 
also leading collective variables even in more complex 
systems. 

When we decrease the density (or saturation) of the LJ 
system, the lifetime (r = 1/fc, k is the transition rate) of 
the gas phase increases drastically. For example, r ~ 10 4 
when the density n = 0.016, but increases by about a 
factor of 20 when the density is only slightly decreased 
to n = 0.014. Even in this density region, direct simu- 
lation of the gas/liquid transition is still possible. Thus, 
we compared the transition kinetics with and without 
the use of a bias. Fig.lJTJ shows that the distribution of 
the gas to liquid transition waiting time t is exponen- 
tial in t from both non-biased and biased simulations, 
lnP(i) cx — t/r. The boost factor a resulting from the 
hyperdynamics method, which characterizes the average 
gain in the rate at which time advances relative to di- 
rect MD, is about 21, as shown in the inset of Fig.JTJ 
(a = 1 for unbiased MD simulation). We also directly 
calculate r by averaging the transition time of full sim- 
ulations for a better comparison. We found that in the 
biased MD case, r = 2.02 x 10 5 , in excellent agreement 
with the direct MD result, r = 2.07 x 10 5 . In our cur- 
rent simulations, we gradually increase the bias potential 
until our desired transition can happen in the usual MD 
steps. Thus we need and can apply larger bias potentials 
in the lower-density LJ systems. At n = 0.012, the phase 
transition is detected while a is of order of 100. At still 
lower density, for example, n = 0.01, it is very difficult 
(if not impossible) to observe the transition using direct 
MD simulations. However, with our method we can still 
easily detect the gas/liquid transition. Fig.© shows the 
results for n = 0.008 and N = 1000. The distribution 
P(52) is flatter and broader in the biased simulation, in- 
dicating that the system visits a larger conformational 
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space. The lower panel of Fig.© shows the reproduced 
free energy profiles from the distribution of the biased 
simulation (a w 10 6 ). It greatly agrees with that from a 
non-biased simulation in the region where the direct MD 
is possible. The inset of Fig.(0) shows the distribution 
of samples in the (V, S2) space. The shown dense region 
(52 > —300) which corresponding to the gas phase is bi- 
ased due to its higher distributed density of samples. The 
lower-density region shown in the inset (S2 < —300) cor- 
responds to the transition region ( the liquid phase which 
52 is far smaller —300 does not shown) was not biased. 
From the obtained distribution, we know 52 is actually a 
good reaction coordinate, the difference between the re- 
built and biased log 10 P(52), shown in the higher panel 
of Fig. (J2J) , corresponds the profile of the applied bias po- 
tential (with a factor fc^TlnlO). 

To summarize, we have expanded the hyperdynamics 
method to more general cases by inhomogeneously com- 
pressing time and conformational space. Our approach 
directly generates an explicit general method to design 
the bias potential. In simple systems, a few functionals 
of the pair correlation function provide a good approx- 
imation of the density field for identifying the impor- 
tant conformations and for constructing the bias poten- 
tial without prior knowledge of the conformational space. 
The method is expected to be applicable in more com- 
plex fluids where even more collective variables might be 
needed. 

This work was supported by the US DOE under con- 
tract W-7405-ENG-36 with LDRD No. LA-UR-06-2794. 
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FIG. 2: (Color online). Top: the distributions of two body 
entropy £2 from non-biased and biased simulations. The re- 
built distribution of the bias simulation is also shown. Here, 
n = 0.008, N = 1000 and T = 0.613. Bottom: the free energy 
profiles from the non-biased and biased simulations are com- 
pared. The inset shows the simulated samples in the (52, V) 
space. The observed liquid phase does not show here. 
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FIG. 1: (Color online). The distribution of transition time of 
the gas phase in direct MD (a = f ) and biased MD (a — 20) 
simulations. The inverse slope of the least-squares fit to the 
points (solid line) gives the lifetime of the gas phase, n = 
0.014, TV = 400 and T = 0.613. Inset: the time boost in a 
typical biased MD simulation is shown. 
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